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Abstract. Starting from delay equations that model field retardation effects, we 
CO ■ study the origin of runaway modes that appear in the solutions of the classical equations 

, of motion involving the radiation reaction force. When retardation effects are small, 

| we argue that the physically significant solutions belong to the so-called slow manifold 

of the system and we identify this invariant manifold with the attractor in the state 
space of the delay equation. We demonstrate via an example that when retardation 
effects arc no longer small, the motion could exhibit bifurcation phenomena that are 
. not contained in the local equations of motion. 
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1. Introduction 

In the treatment of the motion of extended bodies in classical field theory, the 
derivation of radiation reaction forces is based upon certain expansions of the retarded 
field potentials in powers of the retardation The resulting local equations of 

motion involve derivatives of the acceleration and generally suffer from the existence 
of unphysical runaway solutions. Under certain model circumstances, we trace the 
origin of these problems to the expansion of the functions of the retarded arguments 
resulting in the replacement of the original nonlocal delay equations of motion by local 
higher-derivative equations that exhibit runaway solutions. In this general context, the 
properties of the delay equations that appear in classical field theory were first studied 
by L. Bel 0, ^|, |J]. Although our approach is rather general, for the sake of concreteness 
we discuss physical situations involving only the gravitational interaction. 

Consider, for instance, inspiraling compact binaries that are expected to be 
promising sources of gravitational radiation. For a binary that is comprised of two 
compact objects — neutron stars or black holes — with, say, approximately equal masses 
m and m! in nearly circular orbits about each other, the relative orbital radius decays 
because orbital energy is emitted in the form of gravitational radiation. The dynamics 
of a usual binary system can be adequately described using the post-Newtonian 
approximation scheme that is valid in case the gravitational field is everywhere 'weak' 
and the motion is slow, that is v <C c, where v is the characteristic orbital speed and c is 
the speed of light. Although Einstein's equations have a hyperbolic character associated 
with the retarded gravitational interaction, the standard post-Newtonian approximation 
scheme of general relativity deals with functions of instantaneous coordinate time 
t rather than the retarded time t r = t — r/c, where (for the binary system) r is 
the effective distance between the bodies (approximately the relative orbital radius). 
The gravitational potentials, which are originally functions of the retarded time, are 
expanded in Taylor series about t using the effective small parameter in this expansion 
that can be written as co^r/c = v/c, where u;& = 2tt / P& is the relative orbital frequency 
and Pb is the binary period. Because the gravitational waves emitted by the binary 
have an effective frequency of ~ 2u;& and wavelength A& ~ cPb/2, the small parameter in 
the expansion can be reduced to the ratio nr/Xb- Due to the observational fact that in 
typical astronomical systems v/c <C 1, the first few terms of such an expansion can be 
used to derive the post-Newtonian equations of motion that describe, for instance, the 
orbital evolution of the binary pulsars discovered by Hulse and Taylor |5|]. The post- 
Newtonian equations of motion of binary stars are similar to the Abraham-Lorentz form 
in electrodynamics ||, 0, |Sj but, because of the tensorial character of the gravitational 
field, these equations involve not only the third, but the fourth and fifth derivatives 
of the stars' positions with respect to time as well. Schematically, the equation of the 
relative orbital motion reads 

f = F (r) + c~ 2 F 2 (r, r, r) + c" 4 F 4 (r, r, r, r (3) , r (4) ) + c" 5 F 5 (r, r, r, r (3) , r (4) , r (5) ), (1) 
where r is the radius vector connecting the stars, the overdot denotes differentiation with 
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respect to time, r( n ) := d n r/dt n and per reduced mass F (r) = —{G{m + m')/r 3 )r is 
the Newtonian force, F 2 is the post-Newtonian force, F 4 is the post-post-Newtonian 
force and F5 is the gravitational radiation reaction force responsible for the decay 
of the orbital period (Pf, < 0) associated with the emission of gravitational waves 
by the binary system. In the quadrupole approximation under consideration here, 
the gravitational waves carry away energy and angular momentum, but not linear 



momentum || 1C| [11], [12|, [I_3|, |bj, |lj| [U| [17|, [Tj|, |19| ; therefore, the total momentum of 



the binary system is conserved and this fact is responsible for the absence of a force term 
proportional to c -3 in (^. Moreover, all tidal, spin-orbit, and spin-spin interactions are 
neglected in (HD; the only parameters in equation ([!]) are the masses and the separation 
between the centers of mass of the members of the binary system. We mention that 
relativistic hydrodynamical (Euler) equations similar to system (|l[) have been derived 



to describe the motion of the fluid elements of the stars 20, 21 



The higher time-derivative equations of the form ([J) cannot be used directly to 
predict the dynamical evolution of a physical system because of the existence of so-called 
runaway modes that have been much discussed in the literature on electrodynamics 
but not in connection with astrophysical problems involving gravitational radiation 
reaction and the calculation of templates of the gravitational waves emitted by coalescing 
binaries [pO , pi] , p2] , |23| . In analogy with electrodynamics, the existence of these runaway 
modes suggests that the truncated equations of the form (|1]) may not correctly predict 
the qualitative behavior of the solutions of the original true dynamical delay-type 
equations [0] that involve the retarded time t r = t — r/c. Moreover, the existence of 
runaway modes can cause serious difficulties for numerical integration in addition to the 
problems associated with the inaccuracies inherent in the approximation of higher-order 
derivatives by finite differences [EH, EI . 



In case the post-Newtonian expansion parameter r/Aj, is sufficiently small, we will 
provide in the following section a theoretical basis for eliminating the runaway solutions 
by replacing system ([]]) with a new model that is a system of second-order ordinary 
differential equations. Within this theory, high-order vector differential equations like 
equation ([!]) are not the desired approximate equations of motion, and they should 
not be used for numerical integration. Rather, system (H) must be viewed merely as an 
intermediate step in the derivation of the physically correct, second-order model equation 
with no runaway solutions that faithfully approximates the dynamics of the underlying 
delay-type equation. For illustration purposes, we apply this approach in section [3] 
to a discussion of one-dimensional gravitational dynamics of a two-body system. As 
expected, the reduced model predicts the correct motion of the binary except possibly 
when r/Xb is not small and residual terms that have been neglected in equation ([!]) start 
to play a significant role. Indeed, for an inspiraling binary system, the effective delay 
uJbr/c increases to some noticeable finite value as the system approaches coalescence. 
Motivated by this physical scenario, we introduce a simple model involving variable 
delay in section |] that can be expressed as a Duffing-type differential-delay equation. 
This model is then analyzed to show some specific behavior of this delay equation that 
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is not predicted by expansion in powers of the delay. Finally, section [5] contains a 
discussion of our results. 

2. Delay equations with small delays 

The delay-type equations of motion with retarded arguments are usually too complicated 
for mathematical analysis; therefore, we limit our discussion in this section to equations 
with constant delays. Although this is an unrealistic restriction in general, we note that 
for an astrophysical binary system consisting of compact point-like neutron stars or black 
holes moving around each other along circular orbits, the delay is almost constant. In 
fact, a close approximation to this delay is the ratio r/c, where r, the radius of the 
relative orbit, is changing very slowly due to the emission of gravitational energy in the 
form of gravitational waves. 

Taking into account the last remark, let us consider a family of delay differential 
equations of the form 

x(t)=F(x(t-T),x(t)), (2) 

where r is viewed as a real dimensionless parameter and x is a variable in M. n ; intuitively, 
the constant delay r corresponds in effect to u^r/c. The members of this family 
are examples of a more general and widely studied class called retarded functional 
differential equations (see [^4], [2^ 



Using the delay equation (|2|) as an abstraction for the retarded-time model that is 
supposed to be approximated by a system of the form ([!]), we will discuss an approach for 
extracting the 'correct' dynamical equations of motion from system (|I]) that eliminates 
the runaway solutions. 

Our approach assumes the existence of an attractor for the underlying delay-type 
equation. We will rely on the work of Bel || for (numerical) evidence in favor of the 
existence of attractors in the retarded equations of motion with space-dependent delays 
that appear in electrodynamics; but, we know of no mathematical proof for the existence 
of attractors for these equations or for the similar delay-type equations of astrophysics. 
Indeed, the proof of the existence of attractors for delay-type equations with space- 
dependent delays remains a challenging mathematical problem of physical significance. 
For the delay equation (Q), however, if \t\ is sufficiently small, then the corresponding 
member of the family (|2]) has a global n-dimensional attractor such that the restriction 
of the delay equation to this attractor is equivalent to a first-order system Sa of ordinary 
differential equations. We will eventually outline a proof of this result. But, let us first 
discuss our approach to eliminating the runaway solutions. 

The solutions of the delay equation (fj) approach the attractor exponentially 
fast; therefore, the system Sa on this attractor determines (asymptotically) the true 
dynamical behavior of the system, hence we consider it to be the 'correct' physical 
model. On the other hand, it is easy to see that if equation (|2J) is expanded in the small 
parameter r and truncated at some order N, then an iVth-order ordinary differential 
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system S^ akin to system ([!]) is obtained such that the coefficient of the iVth-order time 
derivative of x contains the factor t n and is therefore singular in the limit as r — > 0. 

For t > and sufficiently small, the high-order differential equation Sn has an 
equivalent first-order system S that has an n-dimensional (invariant) slow manifold. 
Moreover, this slow manifold has corresponding stable and unstable manifolds; in effect, 
the first-order system has (physical) solutions that are asymptotically attracted to the 
slow manifold and (unphysical or runaway) solutions that are asymptotically repelled 
from the slow manifold. The restriction of the first-order singularly perturbed system of 
differential equations to its slow manifold is of course an n-dimensional first-order system 
of ordinary differential equations Ss on this n-dimensional manifold. Our main result 
states that in appropriate local coordinates, the system Ss on the slow manifold agrees to 
order N inr with the first-order system Sa on the global attractor of the underlying delay 
differential equation; therefore, the system Ss, which can be obtained directly from the 
high-order differential equation S^, is a faithful approximation of the 'correct' physical 
model. Generalizing to the Abraham-Lorentz type equation ([!]) (analogous to Sn), 
the correct physical model is obtained as the system of ordinary differential equations 
(analogous to Ss) that determines the motion on the slow manifold of a corresponding 
first-order system that is viewed as being singularly perturbed relative to the small 
parameter r/A&. 

While there is evidence that the mathematical assertions in the scenario just 
proposed are valid, some of these assertions have not yet been rigorously justified in 
full generality, even for the case of fixed delays. In the remainder of this section we will 
provide some evidence, in the case of fixed delays, for the existence of a global attractor 
and for the claim that the dynamical system on this attractor is well approximated by 
the dynamical systems on the slow manifolds of singularly perturbed first-order systems 
obtained by truncations of the expansion of the delay equation in powers of the delay. 

Our approach for the elimination of runaway solutions is equivalent to the procedure 
of iterative reduction (also called order reduction) that is often used to eliminate runaway 
solutions by means of the evaluation of the higher time-derivative terms in equations 
like (0) by the repeated substitution of the equations of motion and the subsequent 
reduction of the resulting equation to one of the second order (cf. || ^, |8j). Thus, 
our approach provides a theoretical framework for the rigorous justification of iterative 
reduction (cf. [§, 0, §]), a procedure that has been justified so far by physical intuition. 

For conservative higher time-derivative systems, the order reduction procedure has 
been investigated within the frameworks of Lagrangian and Hamiltonian dynamics by 



a number of authors (see p6|, [27], |28], |29], [3(J an d the references cited therein). In 



particular, it can be shown that under suitable conditions — relevant, for instance, to 
Euler-Lagrange equations analogous to equation ([!]) truncated at the fourth order- 
higher- derivative Lagrangians can be iteratively reduced by redefinitions of position 
variables . 

Returning to the delay equation (^), we note that it has an infinite-dimensional 
state space of initial conditions. For example, if the delay r is a fixed positive number, 
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then the natural state space of initial conditions is the infinite-dimensional vector space 
of continuous IR n - valued functions on the interval [— r, 0]. This space endowed with the 
supremum norm is a Banach space that we will denote by C. Note that for an arbitrary 
continuous R n - valued function 7 defined on the interval [t — r,t], the function j t given 
by 7t(i?) = "fit + "&) is in C. Under the assumption that F is a smooth function and 
G C, there is a unique continuous solution y of the corresponding delay equation in 
the family (§) such that y is uniquely defined for t > — r and y = (see, for example, 
|24|, |25f). The state of the system at time t > is defined to be the function y t in C. 

To see that there is an attractor for the family (0) in case r is sufficiently small, it is 
convenient to introduce the fast time s := t/r, valid for r 7^ 0, so that with y(s) := x(t) 
the family (FJ) takes the form 

y'(s)=rF(y(s-l),y(s)) (3) 

and each member of this family, parametrized by r, has the same state space — the 
continuous functions on the interval [—1,0]. For each r ^ the delay equation fl3|) is 
equivalent to the corresponding member of the family (|2|). For r = the corresponding 
differential equations are not equivalent, but this is of no consequence because we are 
only interested in the solutions of the family (|2]) for r / 0. By viewing the unperturbed 
system (^), namely y'(s) = 0, as a delay equation with unit delay, it is clear that 
the solution with initial state is given by y = on the interval [—1,0] and by the 
constant y = 0(0) for t > 0. The initial state in C thus evolves at time t — 1 to its 
final constant state, the function defined on the interval [—1,0] with the constant value 
0(0). Thus, we conclude that the n-dimensional space of constant functions on [—1,0] 
is a global attractor for the delay equation y'(s) = 0. Moreover, the convergence to 
this attractor is faster than any exponential (the solution reaches the attractor in finite 
time), and the dynamical system on this attractor is given by the ordinary differential 
equation y'(s) = 0. If F is appropriately bounded and r is sufficiently small, then, 
because the contraction rate to the attractor is exponentially fast, the attractor persists 
in the family (^) in analogy with the persistence of attractors in finite-dimensional 
dynamical systems; in fact, each corresponding member of the family (Q) has an n- 
dimensional attractor in the state space C and the restriction of the dynamical system 
to this attractor is an ordinary differential equation. In particular, the family (|3]) has a 
corresponding family of invariant manifolds (that is, manifolds consisting of a union of 
solutions) that depend smoothly on the parameter r. 

To identify the dynamical system on an attractor of a delay equation, let us 
suppose that the delay equation (|2|) has a family of n-dimensional invariant manifolds 
parametrized by r. Moreover, let £ denote the local coordinate on these invariant 
manifolds, and let x(t, £, r) denote the solution with the initial condition x(0, £, t) = £ on 
the invariant manifold corresponding to the parameter value r. Because these solutions 
satisfy the delay equation (0), we have that 

x{t, £, r) = F{x(t - r, £, r),x(t, £, r)); (4) 
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therefore, the generator of the dynamical system on the attractor is the vector field 

X(£,t) := J^(a,r)| t=0 = F(x(-r,e,r),0- (5) 

Under our assumption that this vector field is analytic in r, it can be expanded as a 
Taylor series about r = by differentiating the function F(x(—t,^,t),^) with respect 
to t. To this end, we note that the partial derivatives of x(—t,£,t) with respect to its 
first argument can be evaluated using equation (f|), and partial derivatives with respect 
to its third argument vanish at r = since x(0, £,r) = £; moreover, its mixed partial 
derivatives can be evaluated by differentiation of equation (^) with respect to r. 

As a concrete and instructive example of the construction of the dynamical system 
on an attractor, let us consider a simple case of equation (0) by replacing it with 
x(t) = f(x(t — r)) where the function / : R — > R is the scalar linear function given 
by f( x ) — ax so that the associated family of delay equations is 

x(t) =ax(t-r). (6) 

In this case, there is a corresponding family of solutions given by 

x(t,e,r)=e A ^, (7) 

where A(r) is the unique real root of the equation A = aexp(— At), a fact that is easily 
checked by direct substitution of equation (0) into the delay equation @. 

The dynamical system on the invariant manifold is generated by the family of 
vector fields X(£,r) = f(x(—T,£,r)) = ax{— r, £,r) = aexp(— A(r)r)£ = A(r)£. By 
an application of the Lagrange inversion formula ||31|1 , the Taylor series expansion of 
X(£, r) about r = is 

_°° ( n-l n 

' i (8) 

1 

n=l 

and its radius of convergence is r* = (|a|e) _1 . The qualitative behavior of solutions 
of the system (0) for small r is clear: all solutions are attracted to a one-dimensional 
invariant manifold on which the dynamical system is the linear ordinary differential 
equation x = X(x,t). For example, if a < and r is sufficiently small, then all 
solutions are attracted to the trivial solution x(t) = 0. 

Let us now turn to the standard approach in physics where an underlying delay 
equation is expanded in powers of the delay to obtain an ordinary differential equation of 
motion. To illustrate this, let us consider a special scalar case of the delay equation (§) 
given by 

x = f{x(t-r))+g(x(t)), (9) 

and let us suppose, in analogy with the true dynamical delay-type equations that might 
arise in theories of electromagnetism and gravitation, that the true equation of motion 
for some process is the delay equation (§). The result of expanding equation @ to order 
t 2 is the second-order differential equation (an analogue of equation ([]])) 

x = f(x) + g(x) - T f{x)x + y [f"(x)x 2 + f'(x)x], (10) 
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where a prime denotes differentiation with respect to x. Although we only write the 
second-order expansion, we note that the coefficient of the iVth-order time derivative 
of x in the iVth-order expansion is (r /N\)f'(x). Hence, the corresponding iVth-order 
ordinary differential equation is singular in the limit as r — > 0. Also, if we assume that 
system (|9]) has a (smooth) family of attractors parametrized by r, then the corresponding 
family of vector fields generating the dynamical systems on these attractors is given to 
second order in r by 

X(x, r) = f(x) + g(x) - rf'(x)(f(x) + g(x)) 

+ lr 2 {f'(x)[f(x)+g(x)r 

+ f'(x)[3f'(x) + g'(x)][f(x) + g(x)}} + 0(r 3 ). (11) 

The 'correct' model (that is, the dynamical system on the attractor in the original 
delay equation) can be obtained by treating the expanded and truncated system akin to 



system ([TOD as a singular perturbation problem, which can be analyzed using Fenichel's 
geometric theory of singular perturbations A basic result of this theory states 

that if an iVth-order singular perturbation problem with small parameter r is recast as 
a first-order ('fast') system and the corresponding unperturbed system has an invariant 
manifold that satisfies certain conditions (normal hyperbolicity) , then for sufficiently 
small t each member of the family of perturbed first-order systems has an invariant 
slow manifold. The dynamical system on this slow manifold for the perturbed first- 
order family obtained from the iVth-order truncation of the delay equation is the desired 
faithful approximation to the correct model. For example, let us recast the second-order 
ordinary differential equation ( |i0|) as the first-order singular perturbation problem 

x = u, 

r 2 u = (f(a;))- 1 [2(l + rf'(x))u - 2f(x) - 2g(x) - r 2 f"(x)u 2 ]. (12) 

Using Fenichel's theory, it is easy to show that each member of this family, corresponding 
to a sufficiently small value of |r|, has a slow manifold. Also, it is possible to prove that 
the family of vector fields on these manifolds agrees to order r 2 with the family (p|) 
of vector fields on the attractor in the state space of the underlying family of delay 
equations ®. We note that these results are valid for the vector case of delay 
equation (§) as well. 

For the delay equation (H), and also for more general families of delay equations 
where the delay is viewed as a small parameter, we conjecture that the slow vector 
field, for an appropriately defined first-order system that is equivalent to the iVth-order 
truncation of the expansion of the family in powers of the small delay, agrees to order N 
with the vector field on the attractor in the state space of the original delay equation. 
We have just mentioned that this conjecture is true for the delay equation (|9]) in case 
N = 2. It can be shown that the conjecture is true in general for the linear delay 
equation x(t) = Ax(t — r), where x is a variable in MJ 1 and A is a nonsingular n x n 
matrix GT 
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As we have already discussed, singular equations of motion like system (P^) 
generally have unphysical runaway solutions. To eliminate these solutions and leave only 
the physical solutions, the singular system must be replaced by the dynamical system on 
the corresponding slow manifold. In effect, the truncated equations obtained from the 
underlying delay equation after expansion in the small delay must be replaced by the 
system obtained using iterative reduction; this system is equivalent to the dynamical 
system on the slow manifold. Without this replacement, the appearance of spurious 
runaway modes in inevitable, and their existence will cause overflows in numerical 
simulations. 



3. Gravitational radiation damping 

To illustrate the singular perturbation procedure described in section § as a method for 
the elimination of runaway solutions, we examine an application of this approach to a 
one-dimensional Abraham-Lorentz equation of the form ([]]). 

Let us consider an ideal linear quadrupole oscillator (that is, two masses m and 
wl connected by a spring of negligible mass), where the only source of damping is 
the gravitational radiation reaction force associated with the emission of gravitational 
radiation due to the variable quadrupole moment of the system. A model for the 
(dimensionless) relative position z of these particles, with gravitational radiation 
damping included, is the fifth-order ordinary differential equation 

d 5 z 2 d 2 z , . . 

» z ^r + dP +z = 1 > (13) 

where the small parameter is given by /i = 4G/i ^o^o/(15c 5 ), fio is the reduced mass 
(fj,Q 1 = m _1 + m /_1 ), £ is the spatial scale parameter and u is the frequency of the 
ideal linear oscillator such that Uq 1 is the temporal scale parameter. In equation (|T3|), 
we have neglected the Newtonian gravitational interaction between m and m' as well 
as all relativistic effects except for radiation damping. Let us note that this oscillator 
has an equilibrium solution z(t) = 1. According to our general approach described in 
section |2|, this system can be treated as a singular perturbation problem, where the 
physically correct dynamical system would be the system defined on the slow manifold 
of a corresponding, appropriately chosen, first-order system. 

We emphasize that the fifth-order differential equation (O) is not the correct 
physical model; for example, to specify a solution, the initial position and its first 
four time derivatives must be given. Even with the obvious choice for these initial 
conditions — that is, the initial conditions for a sinusoidal oscillation — a numerical 
integration shows that such solutions do not oscillate; rather, they are divergent. 

To obtain a system with the expected dynamical behavior of an under-damped 
oscillator, the differential equation (|i~3l) must be replaced by its restriction to an 
appropriate slow manifold. We will not carry out the complete reduction procedure 



here fl34|. We note, however, that the system matrix for the linearization of system (|13|) 
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at the steady state solution z = 1 has five distinct eigenvalues that are given to lowest 
order in the small parameter by 

For small /i, the first three eigenvalues are 'fast' and the last two are 'slow'. This suggests 
that the nonlinear system has a two-dimensional slow manifold. In fact, in accordance 
with our general scheme, the restriction of the dynamical system to this invariant 
manifold is a second-order system that gives the correct post-Newtonian dynamics. 
In this case, the dynamical system on the slow manifold to first order in the small 
parameter is given by the second-order differential equation 

15 

z + 32fiz(z )z + z = l. (14) 

16 

A unique solution of this equation is obtained by specifying only the initial relative 
position and velocity of the oscillating masses. For z near the equilibrium state z — 1, the 
expected dynamics for the radiating system is revealed: the relative motion is an under- 
damped oscillator. Numerical integration of this equation using standard algorithms is 
stable and produces the expected result. 

The iterative reduction procedure can also be used to obtain equation ( |T4|) from 
equation ([13]). Even our simple example illustrates the necessity of reducing the higher- 
order equations of motion involving radiation reaction before numerical integration. For 
the more realistic hydrodynamic equations that include conservative post-Newtonian 
terms as well as radiation reaction, the corresponding Euler equation must involve these 
forces in the reduced form, that is, they should contain at most the position and the 
velocity of the fluid element pUL pll . 



4. Delay equations with sufficiently large delays 

It is important to point out that the reduction procedure described in sections |2] and ^| 
cannot in general be expected to produce a good approximation to the true dynamics 
for 'large' delays. 

As a simple but revealing example, let us reconsider the scalar linear delay equation 
x(t) = —ax(t — r) with a > 0. For small |r|, we have already shown that all orbits in the 
state space are attracted to a one- dimensional attractor on which the dynamical system 
is given by the vector field (^|) with a i— > —a. For |r| less than the radius of convergence 
of this series r* = (|a|e) _1 , the correct dynamical behavior of the delay equation is 
predicted by this vector field. Because, in this case, the zeroth-order approximation 
x = —ax already has a hyperbolic structure (that is, all solutions are attracted to 
the rest point at the origin exponentially fast), even the zeroth order approximation 
determines the qualitative dynamics for these values of r. By inspection of this delay 
equation, it might seem natural to conclude that the fixed delay r does not influence 
the behavior for sufficiently large t and the approximation x = —ax remains valid for 
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all fixed delays. This is not true. For instance, if r = ir/(2a), then the delay equation 
has the two-parameter family of exact solutions 

t 1— > c\ cos at + C2 sin at. (15) 

Therefore, the qualitative behavior of the delay equation x(t) = —ax(t — 7r/(2a)) 
is certainly not predicted by the ordinary differential equation x = —ax, or by the 
corrections to this equation within the radius of convergence of the slow vector field. 
The transition of the dynamical behavior of this delay equation from a stable rest 
point to a periodic regime as r increases is easily seen to be the result of a degenerate 



Hopf bifurcation |35]. Indeed, we recall that x(t) = exp(A(r)t) is a solution of the 
delay equation under consideration if A is a solution of the characteristic equation 
A = — aexp(— Ar). For r < r*, the solutions of this equation have negative real parts 
and all such solutions are therefore attracted to the zero solution. If r = 7r/(2a), then 
the characteristic equation has a pair of pure imaginary roots that give rise to the two- 
parameter family of periodic solutions (|T5|). For r > it /(2a), the characteristic equation 
has roots with positive real parts; therefore, there are solutions that grow without 
bound. Nevertheless, for these values of r, the delay equation has an attractor. In fact, 
for 7r/ (2a) < r < 7r/(2a) +27r/a, there is a two-dimensional attractor and the dynamical 
system on the attractor has the form x — 28x + (8 2 + ip 2 )x = corresponding to the 
roots A± = 6(t) ±z</?(t) of the characteristic equation A = aexp(— Ar) with positive real 
parts. As r increases further, the dimension of the attractor increases disco ntinuously 

by two at each r = n/ (2a) + 2iW/a, where N — 2, 3, 4 

The Hopf bifurcation for delay equations with constant delays has been studied in 
detail. For instance, a more sophisticated analysis (see, for example, (§5), p. 341]) shows 
that t = 7r/2 is a supercritical Hopf bifurcation value for the nonlinear scalar delay 
equation 

x(t) = -[l + x(t)]x(t - t). 

Moreover, this system has a nontrivial periodic orbit for each r > 7r/2. 

The delay-type equations of astrophysics generally do not have constant delays. But 
as we have mentioned, for two coalescing neutron stars with nearly equal masses and on 
nearly circular orbits, the delays involved are almost constant; in fact, this 'fast' periodic 
motion evolves as a result of radiation damping on a timescale that is much longer than 
Pb- During this 'slow' evolution, the delay increases as the radius of the binary decreases 
due to the emission of gravitational waves. Motivated by this astrophysical scenario, 
we have studied an oscillator model with a time- dependent delay. This example is not 
intended to be a realistic model, rather it is meant to illustrate some of the bifurcation 
phenomena that occur in delay equations with time-dependent delays. Our example is 
the second-order differential-delay equation 

x{t) + tt 2 x(t) + ax(t -w)- (3x 3 (t - w ) = 0, (16) 

where a, j3 and Q are constant system parameters and x is viewed as the state variable 
of a (Duffing) oscillator with variable delay w(t) such that w{t) + pw(t) = pv. Here p 



Delay Equations 



12 




_j 5 I I I I I I I I I 

100 200 300 400 500 600 700 800 900 

Figure 1. Plot of x versus t for the delay-differential equation (|l6|). Here x(t) = 0.5 
on the interval t < 0, w(0) = 0, a = 0.1,f3 = 0.1, fi = 1, p = 0.006 and v = 10.5. 

and v are constants; hence, w{t) — v is an exponentially decreasing or increasing function 
of time depending on whether p is positive or negative, respectively. In any case, we 
have a dynamic delay that is asymptotic to the constant value v. Note that if p = 0, 
the delay is constant; in this case, the corresponding second-order differential equation 
on the slow manifold (to first order in w) is given by 

x + w(3(3x 2 -a)x+ (a + £l 2 )x - (3x 3 = 0, (17) 

a form of van der Pol's equation. In case Q ^ 0, this differential equation typically has 
a stable limit cycle for w > 0. But for Q = (that is when all forces are retarded), it is 
easy to prove that no periodic orbits exist and most solutions are unbounded. 

A typical plot of x versus t for system ([il]) for Q 2 — a > is given in figure [l], where 
the delay increases from its initial value w(0) — to v — 10.5. The initial response 
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of the system (where the delay is small) is characterized by an oscillation as expected 
from equation (|T7|), which follows from the expansion of equation ([X6|) to first order in 
w C 1. But as w increases, the qualitative behavior of the system is affected by three 
additional bifurcations not accounted for by equation fllT|). At the third bifurcation, the 
stable oscillation disappears. Additional bifurcations of the same type occur if v is set 
to a larger value. Numerical experiments suggest that these bifurcations are not Hopf 
bifurcations; instead, they are 'center bifurcations', where at some parameter value there 
is a rest point of center type and one of the periodic orbits surrounding this rest point 
continues to exist as the parameter is changed. The family (|17]) with the parameter 
values as in figure [l] has a bifurcation of this type as w increases through w = 0. 

The behavior depicted in figure || is suggested by an analysis of the roots of the 
characteristic equation, 

C 2 + n 2 + ae~ Cw = 0, 

for the linearization of the delay equation 
to the existence of centers) are given by 

w e = M^ 2 + acos£vr)- 1/2 , (18) 

where i is a non-negative integer. These are the values of w such that the characteristic 
equation has pure imaginary roots. A computation shows that if i is even, then as w 
increases a pair of pure imaginary roots crosses the imaginary axis into the right half- 
plane, and if I is odd, then the roots cross into the left half-plane. Under the assumption 
that the bifurcations are supercritical, a stable limit cycle appears after the bifurcation 
in the first case; in the second stable limit cycle disappears. For the parameter 

values used to obtain figure |l], the bifurcation values computed from equation are 
(approximately) 0, 3.3, 6.0, 9.9 for £ = 0, 1, 2, 3 such that wg<v. At w = a limit cycle 
appears, at ui\ ~ 3.3 the limit cycle disappears, and so on. Thus, these bifurcations 
account for the appearance and disappearance of oscillations in figure |l|. We note that 
a similar sequence of bifurcations occurs whenever Q 2 — a > 0. On the other hand, if 
Q 2 — a < (for example if Q = 0), then all bifurcation points correspond to roots of the 
characteristic equation crossing into the right half-plane. In this case, the bifurcations 
can be subcritical. Indeed, for Q = 0, numerical simulations indicate that no limit cycle 
appears. As a result, solutions starting near the unstable rest point become unbounded. 

The slow dynamical system, obtained by reduction from a truncation of an 
expansion of a delay equation in powers of the delay, approximates the dynamics on the 
global attractor of the delay equation as long as the delay is sufficiently small; but, as our 
examples show, the ordinary differential equations obtained by expansion, truncation, 
and reduction cannot be used in general to predict the correct dynamical behavior for 
sufficiently large delays. We have mentioned, for example, that the dimension of the 
attractor of a family of delay equations, parametrized by the delay, can increase in 
dimension so that that the corresponding slow vector field is no longer defined on the 
attractor. But this is not the only possible scenario for the appearance of new attractor 



16). The bifurcation points (corresponding 
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dynamics; for example, the attractor could cease to exist or be a manifold for some 
values of the delay. 

The Abraham-Lorentz type equation ([I]) can be used, after reduction to a slow 
manifold, to predict the relative orbital motion of a relativistic binary system in the 
regime where the delay is sufficiently small. The size of the maximum allowed delay 
would have to be computed on a case-by-case basis using the explicit form of the delay 
equation that models the dynamics of a coalescing pair of neutron stars. The results 
of this section show that for sufficiently large delay the attractor does not in general 
correspond to the slow manifold. The question remains whether such a divergence of 
behaviors could ever occur in the case of retarded equations of classical field theory. 
This is an interesting open problem. 



5. Discussion 



It is expected that interferometric gravitational wave detectors that are presently under 
construction will be able to detect signals from massive coalescing binary systems. For 
the analysis of such forthcoming data, it is important to have theoretically predicted 
wave forms ('templates') for the relevant astrophysical processes. To this end, extensive 
computations are necessary that need to take gravitational radiation reaction into 



account GO, pll E2, [23[ . The standard approach leads to higher time-derivative equations 



that involve runaway modes and inevitably produce incorrect results. 

We have determined the source of the difficulty by investigating delay equations, 
which are essentially nonlocal, and the higher time-derivative equations that are 
obtained by truncations of the (post-Newtonian) expansions in powers of the delay. 
For sufficiently small delays, a proper justification is provided for the usual method of 
replacing terms with higher- derivatives by terms with at most first derivatives using 
repeated substitution of the equations of motion ('iterative reduction'). We have shown 
that in the investigation of the solutions of higher-derivative equations that represent 
phenomena involving radiation reaction, it is essential to reduce such equations to the 
corresponding slow manifolds before numerical analysis. Our work suggests, however, 
that unexpected nonlocal phenomena could occur for sufficiently large delays that cannot 
be predicted using the local equations of motion even after iterative reduction. 
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